



# to manipulate data
library(tidyverse)

# to estimate penalized likelihood models
library(brglm2)

# to cluster standard errors
library(lmtest)
library(sandwich)

# to generate results tables
library(stargazer)


# load data
load("full_data_conditionality_covariates_unwaived.RData")

#### SUBSET DATA 
glm_df <- full_df %>%
  select(iso3c,year,short_term,long_term,previous_short_term,previous_long_term,
         imf_program,resource_conditionality,resource_conditionality_dict_tfidf,
         resource_rents,log_gdp_pc_constant,gdp_growth,field_discovery,
         oil_price,crisis,polity2,left_executive,election,any_war,ideal_point_usa,eiti_member,
         nationalization_oc,wb_extractive_project,working_age_pop,t) %>%
  mutate(resource_conditionality = resource_conditionality * 100) %>%
  group_by(iso3c) %>%
  arrange(year, .by_group = T) %>%
  mutate(across(imf_program:working_age_pop, list(lag = ~ lag(.,1)))) %>%
  ungroup() %>%
  filter(year>1979)



# SHORT TERM
brglm1 <- glm(short_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm1_r <- brglm1 %>% 
  coeftest(vcovHC(brglm1, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm2 <- glm(short_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm2_r <- brglm2 %>% 
  coeftest(vcovHC(brglm2, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm3 <- glm(short_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm3_r <- brglm3 %>% 
  coeftest(vcovHC(brglm3, type = 'HC0', cluster = c('iso3c', 'Summary')))


# LONG TERM
brglm4 <- glm(long_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag +
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm4_r <- brglm4 %>% 
  coeftest(vcovHC(brglm4, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm5 <- glm(long_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag +
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm5_r <- brglm5 %>% 
  coeftest(vcovHC(brglm5, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm6 <- glm(long_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit",control=glm.control(maxit=500))

brglm6_r <- brglm6 %>% 
  coeftest(vcovHC(brglm6, type = 'HC0', cluster = c('iso3c', 'Summary')))

# TABLE G.4 (with clustered standard errors)
stargazer(brglm1_r,brglm2_r,brglm3_r,brglm4_r,brglm5_r,brglm6_r, omit = c("iso3c"), #type = "text",
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Resource Rents (% GDP)",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)


# retrieve n, ll etc
stargazer(brglm1,brglm2,brglm3,brglm4,brglm5,brglm6, omit = c("iso3c"), #type = "text",
          dep.var.labels = c("Short-Term Policy","Long-Term Policy"),
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Resource Rents (% GDP)",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)

